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A shadow of a numerical solution to a chaotic system is an exact solution to the equations of 
motion that remains close to the nnmerical solntion for a long time. In a collisionless n-body sys¬ 
tem, we know that particle motion is governed by the global potential rather than by inter-particle 
interactions. As a result, the trajectory of each individual particle in the system is independently 
shadowable. It is thus meaningful to measure the number of particles that have shadowable trajec¬ 
tories as a function of time. We find that the number of shadowable particles decays exponentially 
with time as and that for e £ [~ 0.2,1] (in units of the local mean inter-particle separation 

h), there is an explicit relationship between the decay constant p., the timestep h of the leapfrog 
integrator, the softening e, and the number of particles N in the simulation. Thus, given N and e, 
it is possible to pre-compute the timestep h necessary to acheive a desired fraction of shadowable 
particles after a given length of simulation time. We demonstrate that a large fraction of particles 
remain shadowable over ~100 crossing times even if particles travel up to about 1/3 of the softening 
length per timestep. However, a sharp decrease in the number of shadowable particles occurs if the 
timestep increases to allow particles to travel further than 1/3 the softening length in one timestep, 
or if the softening is decreased below ~ 0.2h. 

PACS numbers: 02.60.-x, 02.70.-c, 05.45.-a, 05.45.Jn, 95.75.P, 05.40.-a, 98.10.-l-z 


Introduction. Numerical simulation of the softened 
gravitational n-body problem is used to gain insight into 
the formation, evolution and structure of gravitational 
systems ranging from galaxies and clusters of galaxies to 
large-scale stucture of the Universe [|^ ^ . Since such sim¬ 
ulations have been used to invalidate theories |^, estab¬ 
lishing their trustworthiness is critical. These simulations 
have several sources of error, including: the use of sev¬ 
eral orders of magnitude too few particles, or discreteness 
noise; the use of approximate force-computation methods 
(the latter two errors are compared in [^); the use of a 
softened potential; the use of finite-timestep numerical 
integration to evolve the system of ordinary differential 
equations; and machine roundoff error. These errors are 
aggravated by the fact that gravitational n-body systems 
are chaotic and display sensitive dependence on initial 
conditions: two solutions with nearby initial conditions 
diverge exponentially away from each other on about a 
crossing timescale [Q, so that any error results in the 
numerical trajectory diverging exponentially away from 
the exact solution with the same initial conditions. The 
phenomenon has been described {eg., 0) as the “expo¬ 
nential magnification of small errors”, implying the pos¬ 
sibility that trajectories of such simulations are the result 
of nothing but magnified noise. 

Fortunately, the purpose of a softened n-body simula¬ 
tion is not to follow the evolution of a particular choice 
of initial conditions, but instead to sample the evolu¬ 
tion of large systems whose initial conditions are drawn 
from a random distribution. As such, we would likely be 
more than satisfied if our numerical solution closely fol¬ 
lowed the evolution of a nearby set of initial conditions. 
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The study of shadowing provides just such a property: a 
shadow of a numerical, or noisy, solution is an exact so¬ 
lution whose initial conditions and subsequent evolution 
remain nearby, in phase-space, to the numerical solution. 
Thus, a numerical solution that has a shadow is essen¬ 
tially an experimental observation of an exact trajectory 
of the mathematical system being modelled. Although 
this observation does not alleviate errors introduced be¬ 
tween the physical system and the mathematical model 
(such as discreteness noise and force softening), it does 
say that the numerical simulation is faithfully solving the 
mathematical model. In it was shown that a sin¬ 
gle particle moving amongst 99 fixed particles is shadow- 
able for several tens of crossing times, and that glitches 
(the point beyond which a shadow cannot be found) tend 
to occur near close encounters. In [|j, we demonstrated 
that if M > 1 particles move in a softened system with 
100 —M fixed particles, then very few particles encounter 
glitches within the first few tens of crossing times. How¬ 
ever, both of these studies used highly accurate integra¬ 
tors to generate the “noisy” trajectories. Although high 
accuracy is commonly used for simulations of unsoftened 
systems, softened simulations most often use the 2nd- 
order symplectic and time-symmetric leapfrog integrator. 

In this paper, we use the leapfrog integrator to gener¬ 
ate noisy trajectories of systems that have M particles 
moving and interacting in a softened potential amongst 
a background oi N — M fixed particles I, I, 0. We 
use normalized units Q in which each particle has mass 
1/Al, and the system diameter, crossing time and average 
velocity all have order unity. We then lead the reader 
through the following observations. First, we observe 
that glitches in the trajectory of a single particle occur 
as a Poisson process (Figure 0). Next, we demonstrate 
that as M increases, shadow durations scale roughly as 
1/M°'®. (The physical significance of 0.8 is unclear and 
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FIG. 1: Histogram of shadow durations of 1000 systems in 
which 1 particle moves in the potential of 99 fixed particles 
with softening e = 0.054 or 1/4 of the mean inter-particle 
spacing. Noisy trajectories were generated using a leapfrog 
integrator with timestep h — 0.011 = e/5. After an initial 
transient, the distribution hts an exponential curve with a 
mean shadow duration of 280 crossing times, indicating that 
the moving particle encounters glitches as a Poisson process 
with a glitch probability of 0.36% per crossing time. 


may be dependent upon other parameters such as N, the 
softening e, and the timestep h.) More importantly, this 
scaling can be experimentally reproduced by superim¬ 
posing trajectories of M single-moving-particle systems 
and taking the shortest shadow of those M systems. In 
other words, particles appear to encounter glitches inde¬ 
pendently of one another (Figure H). Now, the glitching 
and subsequent errant behaviour of just one particle (the 
first to undergo a glitch) in a large simulation is unlikely 
to have a large effect on the reliability of that simulation; 
in fact, as long as only a small fraction of particles have 
glitched, then the reliability of the simulation probably 
remains high. Then, assuming that particles encounter 
glitches independently of one another, we can use the 
distribution of shadow durations of M = 1 systems to 
predict the fraction of shadowable particles as a function 
of time. We find that this fraction is a decwing exponen¬ 
tial with some exponent p. (Figures § and ^). Finally, we 
demonstrate an explicit relationship between fi, N, e, and 
h that holds as long as e is in the range ^ [0.2,1] times 
the mean inter-particle separation and h < e/3 (Figure 
H). This means that given N, e, and the expected dura¬ 
tion of the simulation, one can pre-compute the timestep 
h necessary to have a desired fraction of shadowable par¬ 
ticles remain at the end of the simulation. 

Results. Figure introduces a histogram of shadow 
durations for 1000 softened systems with N = 100, M = 
1. After an initial transient (explained in the discus¬ 
sion of Figure 3]), the distribution fits an exponential 
curve, suggesting that glitches occur as a Poisson process. 

Figure ^ introduces how the average shadow duration 
scales as the number of moving particles is increased. 
For various values of M, we perform 40 experiments in 
which M particles move and interact amongst 100 — M 
fixed pariticles, and plot the mean and standard devia¬ 
tion of the shadow durations. We make the following ob¬ 
servations: (1) a glitch in the local 6-dimensional phase- 
space trajectory of any one particle will cause a glitch 
in the full 6M-dimensional phase-space trajectory of the 
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FIG. 2: How the average shadow duration scales as the num¬ 
ber of moving particles M is increased. Each system is identi¬ 
cal to that described in Figure except now M takes on the 
values 1,3,5,7,10,15,20,25,35. The dots represent sample 
shadow durations, 40 samples for each M, while the “sam¬ 
ple average” is plotted with sample error bars of full width 
la. The “predicted average” is artificially constructed for 
each M = 1,..., 35 W superimposing M samples chosen at 
random from Figure [u and taking the minimum shadow du¬ 
ration of those samples. This demonstrates that the aver¬ 
age shadow duration of an M-moving-particle system is well- 
predicted using M single-moving-particle systems, and sug¬ 
gests that particles encounter glitches independently of one 
another; 280/M°'® is plotted for comparison, although the 
physical significance of 0.8 is unclear. 

M-moving particle system. (2) in a large collisionless 
system, the gravitational potential is governed more W 
the global potential than by inter-particle interactions S , 
and so it is reasonable to expect that particles encounter 
glitches independently of each other. So, perhaps the 
mean shadow duration of an M-moving-particle system 
can be predicted by the mean shadow duration of a sys¬ 
tem with M completely uncoupled 1-moving-particle sys¬ 
tems Q . We test this hypothesis with the “predicted av¬ 
erage” of Figure || by taking M samples at random from 
Figure and taking the shortest shadow duration. We 
see that the predicted curve is well within the error bars 
of the “real” M-moving-particle system. Formally, this 
suggests that the average duration before the occurance of 
the first glitch in the system is statistically independent 
of whether particles interact or not. 

Once one particle in the system encounters a glitch, its 
trajectory after that point is incorrect, and it will pre¬ 
sumably start to “infect” the motion of other particles. 
However, by observation (2) above, we can hope that in 
a large collisionless system, one errant particle will, for 
a time, have negligible effect on the trajectories of other 
particles. In fact, we can guess that, for a time, any 
small fraction of errant particles will have little effect on 
the global behaviour of system. The goal would then be 
to minimize, at reasonable cost, the number of particles 
that have glitched by the end of a given simulation. 

We now take as a working assumption that particles 
encounter glitches independently of one another, and that 
the first small fraction of particles that encounter glitches 
have a negligible effect on the others. That is, we re¬ 
interpret Figure to represent the distribution of shadow 
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FIG. 3: The estimated fraction of particles that would be 
shadowed as a function of time, in a system similar to that 
described in Figure |^. except with all particles moving. This 
is computed simply by taking the cumulative distribution 
function F{t) of Figure where F{t) represents the fraction 
of particles that have undergone a glitch, and then plotting 
1-F{t). 

durations for all the individual particles in a single many- 
moving-particle system. Of course, the figure is likely 
only to be valid for a duration much shorter than 800 
crossing times, as the earlier glitched particles “infect” 
the motion of the remainder, but let us assume it is a 
reasonable approximation for some shorter period. This 
allows us, as a first approximation, to estimate the frac¬ 
tion of glitched particles in a real simulation at a given 
time by computing the fraction of one-moving-particle 
systems that have glitched by that time. Figure 1^ plots 
the opposite—the fraction of non-glitched {i.e., shad¬ 
owed particles as a function of time—derived from Fig¬ 
ure □ by taking its cumulative distribution function F{t) 
and “flipping” it to 1 — F{t). As expected from the Pois¬ 
son process in Figure the fraction of shadowable par¬ 
ticles decays exponentially with a rate corresponding to 
a 0.36% glitch probability per particle per crossing time. 

Figure || is interesting, but is of little use because it 
does not tell us how the shape of this curve varies with 
the total number of particles N, the softening e, or the 
timestep h. However, for reasons we will discuss later, 
we have found that if the timestep h is scaled as 

( 1 ) 

then each of Figures 1,2, and 3 are preserved if e is not 
too small. That is, if N and e are changed in a given 
simulation but the initial conditions are drawn from the 
same distribution, then using Equation to scale the 
timestep will preserve the same degree of simulation reli¬ 
ability from the standpoint of shadowing. Intuitively, the 
scaling of oc e is not surprising and is a commonly used 
timestep criterion. The is more surprising, telling 

us we can increase the timestep as N increases at a fixed 
softening; intuitively, this is because the gravitational po¬ 
tential becomes more smooth with increasing N. 

To demonstrate this scaling, we have performed many 
experiments with various values of h, N, and e. Figure 
^ summarizes the results. Each line represents the frac¬ 
tion of shadowable particles as a function of time for some 



FIG. 4: Curves demonstrating that (i) Equation (|^) preserves 
simulation reliability from a shadowing standpoint, (ii) reli¬ 
ability increases as h decreases, and (iii) the scaling breaks 
down if e is too small. Each line represents a curve similar 
to that in Figure ^ derived from a 200-sample set of simula¬ 
tions similar to Figure but with different N and e. There 
are 6 clusters of 4 lines each. The 4 lines in each cluster 
come from 4 sets of simulations using all 4 combinations of 
N = 100,1000 and e = Each system was in¬ 

tegrated using leapfrog with timestep h scaled according to 
Equation (|^. Each cluster represents a particular choice of 
constant k in the scaling, namely, for the 

displayed values of k. Decreasing k increases accuracy, giving 
longer shadows. Embedded in each cluster is a dotted line rep¬ 
resenting discussed later. Finally, the diamond and 

curves have £ = {|, respectively, and demonstrate 

that shadows are much shorter if the softening is too small, 
even though they use the smallest timestep factor k — 


{N, e) pair. Each closely clustered set of 4 lines represents 
sets of runs with various (A^, e) and the the timestep h 
scaled using Equation to give the same “shadowing reli¬ 
ability”. Finally, decreasing the timestep (via decreasing 
k) increases reliability by decreasing the decay rate of 
the fraction of shadowed particles. We also performed 
similar experiments with all combinations of parame¬ 
ters N = 10^,10^,10'* and £ = {2,1, i, i, i, 

The scaling with N works well up to = 10'*, and 
presumably beyond. All curves are very similar for 
e > However, as seen in the figure, shadows 

are significantly shorter for e = even if 

the smallest timestra is used. This is consistent with 
the observations of Q , and may be related to unphysical 
results obtained with a too-small softening 0 - 

Finally, we come to the crux. Armed with the scaling 
Equation and the knowledge that the fraction of shad¬ 
owable particles decays exponentially as we would 

like to find a relationship between the timestep propor¬ 
tionality constant k of Figure and p,k. An eyeball fit 
of exponential curves to each of the clusters in Figure 
^ gives values of the decay constant p,k for each cluster. 
These values, along with the curve /i^ vs. k, is plotted 
in Figure |[ As can be seen, there is evidence that for 
fc < ^, the curve settles to a power law of approximately 
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FIG. 5: For fc = 1, i, i, 4, i, | and e > of Figure Q 

the values of fik were fit by eyeball, and are, respectively, /i = 
0.33, 0.024, 0.0071, 0.0040, 0.0027, 0.0013. These are plotted, 
and for k < fit the curve jj, — (0.047 ± . 


= (0.047 ± 0.015)/c1-^^^°-25. (2) 

However, the shape of the curve is also consistent with a 
slowly decreasing slope as k decreases. 

Discussion. The relations described above offer an 
a priori algorithm for choosing a timestep for a softened 
A^-body simulation, viz. : given N, e, the expected simu¬ 
lation duration T in crossing times, and a desired fraction 
F of shadowable particles remaining at time T, solve for /i 
in F = and then solve for k using Equation (|^). Of 

course, these relationships will need to be scaled to appro¬ 
priate units for the simulation. The job should be easiest 
for a simulation of one galaxy; for simulations of clusters 
of galaxies or a cosmological simulation, we would scale 
to the smallest sub-systems we expect to accurately inte¬ 
grate. We are unsure of the effects of dynamically chang¬ 
ing the softening based on the local mean particle density, 
but suspect that some reasonable interpretation may be 
possible whereby a softening and timestep (modulo the 
discussion of the next paragraph) are chosen based upon 
local mean particle density. 

The fact that constant-timestep leapfrog is symplec- 
tic is probably significant to these results. We exper¬ 


imented briefly with a dynamically changing timestep, 
but found that shadows were virtually destroyed if the 
timestep changes “too often”. However, we found that if 
the timestep was decreased as a particle entered a high- 
density region, but never increased the timestep again, 
these results were preserved. This may be a reasonable 
choice for simulations of clusters if most particles that 
enter a high-density region remain there for the remain¬ 
der of the simulation. Alternatively, perhaps a particle’s 
timestep could be re-increased only after the dynamic 
timestep criterion says the particle’s timestep should be 
significatly increased, say by an order of magnitude. This 
will ensure the particle has left the high-density region 
far behind, and preserve the internal reliability of high- 
density regions. 

More detailed arguments deriving Equation ( 0 ) show 
that the forward global error of a softened Wbody sim¬ 
ulation scales as . The h? scaling is due to 

leapfrog being a globally 2nd-order integrator; the 
scaling has been seen before Q, and the is new. 
Thus, Equation (|^ simply holds the forward global error 
constant. Finally, a shadowing concept known as brittle¬ 
ness relates the forward global error to the distance 
between a shadow and the corresponding points on the 
numerical trajectory, although to our knowledge this pa¬ 
per is the first to demonstrate a relationship between the 
the forward global error and the shadowing distribution. 

Since the scaling of h with is so weak, and we 

usually increase N in order to increase reliability, it cer¬ 
tainly will not hurt to ignore the scaling with N and 
simply scale h (x e. The remaing questions are: what 
value of e to use, and what fraction of e should a particle 
be allowed to travel in one timestep? The first question 
is answered by the fact that shadows do not appear to 
exist for very long if e < concerning the second 

question, we believe that Figures 0 and || suggest that 
the particle should be allowed to travel at most in 
one timestep, and possibly much less, if a shadow dura¬ 
tion of 100 crossing times encompassing most particles is 
desired. 
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